source(paste0(workdir, "scripts/conley_se_files/ConleySEs_17June2015.R"))
# http://darinchristensen.github.io/2015/08/30/conley/
# specify latitude and longitude variables as cluster to have them in output:

d1 %>% filter(smpl == 1) %>% mutate(cnamef2 = cnamef) -> dXc

bf %>% gsub("cnamef$", "caplat + caplong", .) -> xf
runmodel(xf, data = dXc, kcx = T) -> r

haccilow <- haccihi <- hacse <- vector("list", length(r))

for(i in 1:length(r)) { # i = 1
SE <- ConleySEs(reg = r[[i]],
    unit = "cnamef",
    time = "periodf",
    lat = "caplat", lon = "caplong",
    dist_fn = "SH", dist_cutoff = 1000,
    lag_cutoff = 5,
    verbose = FALSE)

# hacse[[i]] <- if (!is.null(nrow(sapply(SE, function(x) diag(sqrt(x)))))) {
hacse[[i]] <- if (nrow(SE[[1]]) == 1) {
    sapply(SE, function(x) diag(sqrt(x)))[3]
} else { sapply(SE, function(x) diag(sqrt(x)))[, 3] }

haccilow[[i]] <- as.numeric(r[[i]]$coefficients - (1.96 * hacse[[i]]))
haccihi[[i]] <- as.numeric(r[[i]]$coefficients + (1.96 * hacse[[i]]))
}

save(r, haccilow, haccihi, file = paste0(workdir, "data/conleyin_d1.RData"))
